// Bibliotecas
#include <cstdio>
#include <iostream>
#include <cstdlib>
#include <string>
#include <fstream>
#include <sstream>
#include <vector>
#include <cmath>
#include <sys/time.h>

using namespace std;

typedef struct vector<double> vect; // vetor de double
typedef struct vector<vect> matrix; // matriz de double
 
// Soma as linhas da Matriz para verificacao de convergecia
double sum_line(matrix& NA, int J_ORDER, int i) {
	double out = .0;

	for(int j = 0; j < J_ORDER; j++) {
		out += NA[i][j];
	}
	
	return out;
}

// Soma do passo de transicao
double sum(matrix& NA, vect& X, int J_ORDER, int i) {
	double out = .0;

	for(int j = 0; j < J_ORDER; j++) {
		if(i != j)  out += NA[i][j] * X[j];
	}
	return out;
}

// Soma a linha do sistema linear para um teste
double sum_ls(matrix& A, vect X, int J_ORDER, int J_ROW_TEST) {
	double out = .0;

	for(int i = 0; i < J_ORDER; i++) {
		out += A[J_ROW_TEST][i]*X[i];
	}

	return out;
}

// 
double max(vect& V, int J_ORDER) {
	double out = 0;

	for(int i = 0; i < J_ORDER; i++) {
		out = out < V[i] ? V[i] : out;
	}

	return out;
}

double erro(vect X, vect& diff, int J_ORDER) {
	for(int i = 0; i < J_ORDER; i++) X[i] = abs(X[i]);
	return max(diff, J_ORDER)/max(X, J_ORDER);
}

double test_convergencia(matrix& NA, int J_ORDER) {
	double max = 0;
	for(int i = 0; i < J_ORDER; i++) {
		double value = sum_line(NA, J_ORDER, i);
		max = max < value ? value : max;
	}
	return max;
}


double jacobi_richardison(matrix& A, vect& B, int J_ORDER, int J_ROW_TEST, double J_ERROR, int J_ITE_MAX) {

	matrix NA;
	vect NB;

	int i, j, k;

	// Matriz *A e *B modificadas
	NA.resize(J_ORDER);
	NB.resize(J_ORDER);
	for(i = 0; i < J_ORDER; i++) {
		NA[i].resize(J_ORDER);
		for(j = 0; j < J_ORDER; j++) {
			if(i == j) NA[i][j] = .0;
			else NA[i][j] = A[i][j]/A[i][i];
		}
		NB[i] = B[i]/A[i][i];
	}

    // teste de convergencia
	if(test_convergencia(NA, J_ORDER) >= 1) return -1;  // caso nao convirga, cai fora
	
	// vector para interacoes
	vect X(NB);
	for(k = 0; k < J_ITE_MAX; k++) {
		vect S(J_ORDER);
		vect diff(J_ORDER);
		for(i = 0; i < J_ORDER; i++) {
			S[i] = NB[i] - sum(NA, X, J_ORDER, i);
			diff[i] = abs( S[i] - X[i] );
		}
		if(erro(X, diff, J_ORDER) <= J_ERROR) break; // Atingiu o erro
		X.clear();
		X = S;
	}
	
	// print results 
	cout << "Iterations: " << k << endl;
	cout << "ROW_TEST: " << J_ROW_TEST << " => "; 
	cout << "[" << sum_ls(A, X, J_ORDER, J_ROW_TEST) << "] aprox " << B[J_ROW_TEST] << endl; 	

	X.clear();

	return k;
}


int main(int argc, char* argv[]) {
	static int J_ORDER;
	static int J_ROW_TEST;
	static double J_ERROR;
	static int J_ITE_MAX;

	// interadores
	int i, j;

	// Matrizes do sistema linear
	matrix A;
	vect B;

	// variavel para pegar as linhas do arquivo que esta em fin
    string line;
    ifstream fin(argv[1]); // ifstream de arquivo
	
	// Leitura do arquivo
	if(fin.is_open()) {
 		string buffer;
        getline(fin, line, '\n'); // pega a linha 1 --> J_ORDER
        buffer = line;

        getline(fin, line, '\n'); // pega a linha 2 --> J_ROW_TEST
        buffer += " " + line;

        getline(fin, line, '\n'); // pega a linha 3 --> J_ERROR
        buffer += " " + line;

        getline(fin, line, '\n'); // pega a linha 4 --> J_ITE_MAX
        buffer += " " + line;	

        istringstream iss(buffer);

       	//DESCARREGANDO O BUFFER
		iss >> J_ORDER;
		iss >> J_ROW_TEST;
		iss >> J_ERROR;
		iss >> J_ITE_MAX;

		// Matriz A
		A.resize(J_ORDER);
        for(i = 0; i < J_ORDER; i++) {
        	A[i].resize(J_ORDER);
        	getline(fin, line, '\n');
    		istringstream iss_line(line);

    		for(j = 0; j < J_ORDER; j++) {
    			iss_line >> A[i][j];
    		}
        } 

        // Matriz B
        B.resize(J_ORDER);
        buffer.clear();
        for(j = 0; j < J_ORDER; j++) {
			getline(fin, line, '\n');
			buffer += " " + line;		
		}
		
		istringstream iss_line(buffer);
        for(j = 0; j < J_ORDER; j++) {
			iss_line >> B[j];
		}

		// Chamada do metodo numerico
		double it = jacobi_richardison(A, B, J_ORDER, J_ROW_TEST, J_ERROR, J_ITE_MAX);
        
		if(it == -1) cout << "not resoveld linear system" << endl;

        // fecha arquivo
        fin.close();

    } else cout << "not found " << argv[1] << endl; // mensagem de erro do arquivo
    
	
	return 0;
} 
